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ABSTRACT 

In this paper we present numerical solutions to several optimal control problems for an 
unsteady viscous flow. The main thrust of this work is devoted to simulation and control 
of an unsteady flow generated by a circular cylinder undergoing rotary motion. By treating 
the rotation rate as a control variable we formulate two optimal control problems and use 
a central difference/pseudospectral transform method to numerically compute the optimal 
control rates. Several types of rotations are considered as potential controls and we show that 
a proper synchronization of forcing frequency with the natural vortex shedding frequency 
can greatly influence the flow. The results here indicate that using moving boundary controls 
for such systems may provide a feasible mechanism for flow control. 
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1. INTRODUCTION 


The benefits received from applying control mechanisms in viscous flows has been long 
realized since the pioneering work of Prandtl [30]. Because of a growing interest in controlling 
the behavior and structure of fluid flows, various topics in flow control have recently become 
a subject of research focus [14, 1, 13, 34, 4, 15, 26, 35, 9, 16, 20, 27, 28, 29, 36]. A potentially 
important application of flow control is the enhancement of the aerodynamic characteristics 
of future advanced aircraft. Although recent research has produced a variety of methods 
to achieve desired performance and design goals for maneuvering vehicles, it is now widely 
recognized that further gains will most likely come from the application of various types 
of flow control mechanisms [11, 12]. Consequently, the utility of flow control becomes a 
critical issue in the design process which may provide real-time effect for many important 
applications, such as highly instantaneous maneuvers for the super-maneuverable aircraft 
[18], and the optimum design of aerodynamic configurations [21]. Considerable effort has 
been devoted to the improvement of control mechanisms. However, the principal progress to- 
date has been essentially accomplished by experimental investigations, while most analytical 
and numerical approaches have remained in its infancy due to the complexity of the problems. 

One of the most practical applications of control mechanisms in flow systems has been 
boundary-layer separation control. In this area of research, several methods have been de- 
veloped experimentally to provide various control mechanisms, for example moving surfaces, 
blowing, suction, injection of a different gas, etc [31, 32]. In particular, it has been demon- 
strated in a number of experiments by Modi et al. [25, 24] that moving surfaces can effectively 
provide boundary-layer control. In their experiments, the boundary-layer flow is controlled 
by an application of two rotating cylinders located at the leading and trailing edges of an 
airfoil. It has been shown that this mechanism can prevent flow separation by retarding the 
initial growth of the boundary layer, with the important consequences of lift enhancement 
and stall delay. For instance, when the speed ratios (which represents the ratio of cylinder 
surface speed to the freestream speed) of both cylinders were set at a constant value of 4, the 
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results indicated a 200% increase of the maximum lift coefficient compared with the reference 
airfoil (in which no rotating cylinder is attached). In spite of the fact that considerable aero- 
dynamic benefits were gained by changing the cylinder speed ratio, in their experiments the 
speed of rotation was performed merely with constant values. However, it should be noted 
that if the rotating cylinder mechanism is applied to a region of unsteady flow, a constant 
rotation rate may not correspond to the optimal performance when an airfoil is undergoing a 
rapid maneuver. This type of result provided the motivation for us to consider a fundamen- 
tal problem regarding unsteady flow control by means of a time-dependent moving surface 
mechanism. In order to keep the problem reasonable and yet practical, we selected a model 
for the numerical study of controlling the temporal development of the flow field around a 
rotating cylinder. 

The most distinguishing feature of a rotating body traveling through a fluid is that it 
experiences a transverse force acting in a direction perpendicular to that of flowing stream 
[37]. In the past few decades, research on the problem of a uniform stream past a cylindrical 
rotating body has been the subject of many experimental and numerical investigations. See 
the papers by Taneda [39], Mo [23] and Tokumaru and Dimotakis [40] for a cylinder under- 
going rotary oscillations, Prandtl [31], Taneda [38], Koromilasand Telionis [22], Coutanceau 
and Menard [8], Badr and Dennis [3], Badr et al. [2] and Chen et al. [7] for a cylinder with 
a constant speed of rotation. However, most of these results are primarily focused on the 
study of formation and development of vortices in a cylinder wake. It appears that the effect 
of the rotation rate on the cylinder forces exerted by the fluid has received far less attention, 
despite the fact that it has many important practical engineering applications. 

The main thrust of the current investigation is on simulation and control of an unsteady 
flow generated by a circular cylinder undergoing a combined (steady or unsteady) rotary 
and rectilinear motion. By treating the rotation rate as a control variable in this model, 
we consider several problems concerning the temporal development of forces on a rotating 
cylinder in response to a variety of time-dependent rotation rates. The computational results 
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provide considerable insight into the problem of controlling such flows. 

2. MATHEMATICAL AND NUMERICAL FORMULATIONS 

In this section the governing equations and the particular numerical method used in this 
work are described. We consider control problems for a two-dimensional viscous incompress- 
ible flow generated by an impulsively started circular cylinder. The cylinder is translated 
with a constant rectilinear speed U normal to its generator and is simultaneously rotated 
with a time-dependent angular velocity Q(t) about its axis. Although there are various 
formulations and numerical techniques for the solving of steady and unsteady flow past a 
rotating cylinder [19, 3, 23], in this work the problem is investigated numerically by solving 
a velocity /vorticity formulation of the Navier-Stokes equations with an implementation of 
the Biot-Savart law. The numerical approach used in the present study is the one developed 
by Chen [5] for the problem of a circular cylinder oscillating in a rectangular box. It is based 
on an explicit finite-difference/pseudo-spectral technique to yield time accurate solutions 
to the governing equations. This numerical algorithm was further modified to investigate 
an unsteady flow around a rotating cylinder undergoing various constant rotational speeds 
[7, 28], and time-dependent rotation rates [26, 27]. 

Several prominent features in the application of an integral representation for flow kine- 
matics are emphasized in §2.2. This integral method proposed by Wu and Thompson [46] 
provides the basic link between the velocity and vorticity fields throughout the numerical 
procedure. The boundary vorticity at the solid surface can be easily calculated by the ap- 
plication of this integral. Moreover, unlike other numerical approaches, the imposition of 
the artificial far-field boundary condition for the velocity is not necessary once the vorticity 
values are known everywhere in the domain of interest. 

2.1 Governing Equations 

In the velocity /vorticity formulation of the Navier-Stokes equations, the governing equa- 
tions consist of the vorticity transport equation and the vector Poisson equation for the veloc- 
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ity. Thus, for a two-dimensional unsteady viscous flow in incompressible fluid, the Cartesian 
coordinate form of the governing equations for the vorticity and velocity can be expressed 
in the dimensionless form as 


du _ 2 2 

- + „.Vo,= -V W 

(1) 

V 2 u = —V x (u>e z ), 

(2) 


where u is the velocity field, u> is the vorticity field, and e z is the unit vector in the direction 
of ^-direction. All the variables are made dimensionless by means of the characteristic 
quantities. The cylinder radius a is used as the length scale while a/U is used as the time 
scale. The Reynolds number Re = 2 Ua/v is based on the cylinder diameter 2a and the 
magnitude U of the rectilinear velocity. 

A non-rotating reference frame, translating with the cylinder is employed. In this frame 
the dimensionless boundary conditions for the problem of a rotating cylinder (with a time- 
dependent angular velocity f l(t)e z ) can be written as 

u = —a(t)ye x + a(t) xe y for ( x,y ) G T (3) 

and 

u = e x for \Jx 2 + y 2 — ► oo, (4) 

where T denotes the impermeable solid boundary of the cylinder. The angular/rectilinear 
speed ratio a(t) = fl(t)a/U is the primary control parameter throughout this work. 

In many practical numerical simulations for the laminar motion of viscous incompressible 
fluid, both the exterior and interior flow problems, the formulation based on velocity /vorticity 
variables would provide some advantages over the primitive- variable formulation. This ve- 
locity /vorticity formulation is especially well suited to treating initial development of the 
flow generated by impulsively started bodies, in which a relatively small vortical viscous re- 
gion is embedded in a much larger inviscid potential flow. Consequently, the computational 
domain may be restricted to a smaller region where all vorticity contributions are contained. 
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Furthermore, the decoupling of the overall problem into its kinematic and kinetic aspects of 
the flow field gives an additional convenience. Throughout this work the vorticity transport 
equation described above may be viewed as kinetic process of the flow field in which the 
distribution of vorticity is interplayed by the process of convection, as well as the effects of 
viscous diffusion. 

2.2 Integral Representation for Flow Kinematics 

Having defined the governing equations for the rotating cylinder problem described 
above, we can now examine several advantages of an integral formulation used in this work. 
In the formulation of exterior flow problems, it is well known that one of the main numerical 
difficulties is related to the proper imposition of a prescribed limit value at the far-field for 
the unbounded physical domain in which the flow takes place. It should be pointed out that 
in any numerical simulation one has to restrict the exterior infinite domain to be finite with 
an artificial boundary. However, the far-field boundary in (4) fails to represent the exact 
characteristics at the outer perimeter of the finite computational domain. Therefore, in many 
practical applications, instead of directly applying the far-field boundary condition (4) one 
often tries to avoid the difficulty by utilizing various asymptotic boundary conditions at large 
distance (e.g. the application of potential flow or Ossen expansion), while others introduce 
a mapping of the infinite domain onto the finite one by means of a suitable coordinate 
transformation. 

The other difficulty encountered in the simulation of viscous flow is that of prescrib- 
ing the appropriate non- velocity boundary conditions at the solid surface. In general, the 
prescribed pressure at the solid boundary is needed in the application of primitive-variable 
(pressure/velocity) formulation, while the boundary vorticity is required for the formulation 
based on the velocity /vorticity (or stream-function/ vorticity) variables. In order to over- 
come these two difficulties, we pose the kinematic relation on the problem by introducing a 
general Biot-Savart induced law described below. 
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In the problem with a viscous fluid, if the velocity distribution of u are given, then the 
vorticity field u could be evaluated through the kinematic relation between u and 2 described 
by 2 = V xu and V • u = 0. On the other hand, the vector Poisson equation 

V 2 u = -V xw, (5) 


obtained from the continuity equation and the definition of vorticity, can be used to determine 
the velocity field from a given vorticity field. 

For a general viscous flow in a region D bounded by an inner boundary T and outer 
boundary T', we can recast the kinematic part of the problem into an equivalent integral 
formulation 


^ 1 f f 2(r, t) x (r - ro) 4 

«(<•„,<) = -JJ d JA 


c J r+r' 


rol 4 

[ub{r,t) • n](r — r 0 ) - [u 6 (f,t) x n] x (r - f 0 ) 
|r-ro|^ 


(6) 


where 


f 4x, ford=3 
\ 2tt, ford = 2. 

In the above integral representation, uj, is the boundary velocity, n is the outward normal 
unit vector and d is the spatial dimension. The subscript “0” denotes the field point where 
the velocity field is evaluated. 

By applying the no-penetration and no-slip conditions to the rotating cylinder problem 
considered in this study, in two dimensions equation (6) can be written in terms of the 
rectilinear velocity Ue x and the angular velocity Q(t)e* z of the solid body B which is known 
as the generalized Biot-Savart law of induced velocity: 

Iff 2(r,t) x (f- ro) 


= ~tJL 


|r — ro| 2 

1_ f f 2 Q (r,t) x (r- r 0 ) 

27T J Jb 


dA 

dA + U. 


( 8 ) 


|r-ro| 2 

For detailed discussions of the integral representation for viscous flows the reader is referred 
to [46, 45]. Equation (8) represents the kinematic relationship between the velocity and vor- 
ticity fields of an infinite domain which is expressed as an integral form. The first integral 
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represents the contribution of the vorticity field to the development of the velocity field over 
the region D occupied by the viscous fluid, while the second integral gives the contribu- 
tion from the rotation of the solid body. Notice that both the no-penetration, the no-slip 
boundary conditions on the solid surface and the far-field boundary conditions are implicitly 
imposed within the integral (8). To be precise, if we let the calculated field point f 0 go to 
infinite in (8), the vorticity contribution will approach zero. Consequently, the boundary 
condition at infinity is satisfied exactly in this integral representation of the kinematic rela- 
tion. This indicates that the difficulty resulting from the imposition of far-field condition is 
removed by an application of (8). 

In numerical simulations, it has been reported that the imposed boundary conditions at a 
large distance from the body will significantly influence the accuracy of the overall numerical 
solution [10, 19]. It is therefore necessary to devise a technique to impose the appropriate 
condition at the outer boundary of the computational domain. The integral approach in 
(8) provides a useful method to accomplish this. To be precise, the integral representation 
(8) permits us to determine the velocity (point-by-point) explicitly if all vorticity values are 
known everywhere in the domain of interest. It will exhibit a more realistic behavior at the 
outer perimeter of computational domain than those asymptotic techniques employed by 
other formulations. Namely, if the computational domain is large enough to contain all of 
vorticity generated around the solid boundary prior to a certain time, then at this instant 
the velocity on the outer perimeter of computational domain can be evaluated directly by 
the numerical integration of (8) with all the known vorticity field in the domain. 

Additionally, if we apply equation (8) to the points of r b on the solid boundary, then the 
integral formula becomes 


1 rr2?i(r,t)x(r-r b ) 

2 nJ Jb |f'-rl| 2 


dA + U. 


The boundary vorticity values are contained in this integral. Hence, by using the prescribed 
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motions of body (i.e. U and Sl(<)) and all of the known vorticity field in the domain, the 
determination of vorticity values at grid points on the solid boundary can be achieved by the 
numerical integration of equation (9). However, many implementations of flow, simulations 
make use of approximate formulae for the boundary vorticity, which may lead to excessive 
error and thereby destroying the accuracy of the solutions. 

When we solve the velocity field in an infinite domain we need only to take the viscous 
region (uj ^ 0) into account, while contributions from the inviscid region vanishes (w = 0). 
Hence, an application of this integral representation allows one to confine the computation 
to the viscous region corresponding to the non-negligible vorticity portion of the flow, in 
contrast to employing domain methods (i.e. the finite-difference and finite-element methods) 
in which both potential and viscous regions are needed in a simultaneous solution. This 
implies the computational domain can be significantly reduced and hence a smaller number 
of grid points are needed than those required by standard finite-difference and finite-element 
methods. 


2.3 General Body-Fitted Coordinate 


In order to accommodate problems with a time-dependent domain in the physical (x, y) 
space, the vorticity transport equation (1) and Poisson equation (2) can be recast in terms 
of the time- varying generalized body-fitted coordinate system (f , rj) as 

w* = J ~ ~ j(“S x v ~ Ur, x () 

-4 [Vr + x v (vu)t - 

+ J^J2 ( au u - ^P^fr + + Q U v)' ( 10 ) 


and 

{ crux - 2f3v.fr + 7 u nr) + J 2 (Pu { -I- Qur,) = J(x n u € - X{U n ), 
avfr - 2 (3vfr + + J 2 {Pv i + Qv n f ) = J(y^ - y^), 

where 

f a = x 2 n + yjj, f3 = x^Xr, + y 7 = x} + yf, 

\ p = + lyy, Q = r] xx + T)yy, J = x^y n - X^Vfr 


(11) 

(12) 
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Here J is the Jacobian of the coordinate mapping between the physical space (x, y) and the 
computational space (£, 77). All subscripts in (10)-(12) denote partial differentiation. Also, 
u and v denote as the velocity component in x- and y-direction, respectively. Notice that 
this coordinate transformation will introduce two additional time-dependent terms x t and 
y t when a time-dependent flow domain in the physical space is considered. However, one 
can calculate the solution in a fixed time-independent computational grid in this generalized 
coordinate system because of the inclusion of x t and y t in the governing equations. Therefore, 
even when a time-dependent domain is considered, the interpolation of boundary conditions 
and grid points employed in the physical space formulation is not necessary in this generalized 
coordinate system. Moreover, it is also convenient to distribute the grid points within the 
domain of interest if a complex flow geometry is encountered. In fact, such a distribution of 
grid points is particularly important for simulations in viscous flows, where the grid points 
clustering near the surface is necessary to resolve the large gradients which appeared in a 
thin region due to the viscous effect. 

Because we wish to extend the results presented here to more general problems with 
time-dependent domains, the computer code was written under this general body-fitted 
coordinate formulation. However, in this work, the polar coordinate (r, 9 ) is used for the 
study of an unsteady flow around a rotating cylinder and the flow domain considered here 
is time-independent since a reference frame translating with the cylinder is used. If a body 
of arbitrary shape other than circular cylinder is considered (e.g. a rotating elliptic cylinder 
or airfoil), then the domain of interest becomes a time-dependent feature even though the 
reference frame is translating with the body. This occurs because of the rotational motion 
of an asymmetric bodies with respect to a non-rotating reference frame. 

2.4 Computational Procedure 

The computational domain is discretized in the body-fitted coordinate system (£, 77) by 
setting & = z'A£, rjj = jAr 7, and i = j = 1 where M and N denote the 

number of grid points in the £- and 77-direction, respectively. The numerical approximation of 
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the vorticity at the grid point (&, T}f) and time t n = nAt is denoted by = u n (iA£,jAr]). 
The vorticity transport equation (10) is first discretized by a second order central difference in 
the radial direction and a pseudospectral transform method in the circumferential direction 
for all spatial derivatives. This semi-discretization form of equation (10), consisting of a 
system of ordinary differential equations in time can be written as 

— = F(u>), u = (u 2 ,2,-” (13) 


for all the interior grid points. Therefore, the calculation procedure to advance the solution 
for any given time increment can be summarized as follows: 

Step 1: Internal vorticity over the fluid region at each interior field point is calculated by 
solving the discretized vorticity transport equation. An explicit second-order rational Runge- 
Kutta marching scheme based on the work of [41] is used to advance in time for (13) 

The discretization of (13) in time thus can be written as 


£ n+1 = Cj n + 


2^1(91,93) ~ 93(91,91) 

(93,93) 


(14) 


( g\ = F(u n )At 

< g 2 = F(u n +cg 1 )At (15) 

l 93 = (1 - h)gi ~ bg 2 

where ( : , : ) denotes the scalar product. In order to ensure the stability of the above nonlinear 
explicit scheme, the two constants b and c in (15) must satisfy be = —0.5. In particular, 
b = — 1 and c = 0.5 are used in our computations. Although this method is explicit in 
nature, it may become unconditionally stable by the suitable choice of the constants in (15) 
[17]. In addition, this particular scheme allows one to use a larger time-step than that of the 
three-step Adams-Bashforth scheme used by Chen [5]. This step consists of the kinetic part 
of the computational loop. 

Step 2: Using known internal vorticity values at all the interior grid points from step 
1, the generalized Biot-Savart law of induced velocity (9) is used to update the boundary 
vorticity values at all the surface nodes. 
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Step 3: At this stage, all the vorticity values in the computational domain are known at 
the new time level. Then, the velocity at points on the outer perimeter of the computational 
domain is calculated by the integral kinematic constraint (8). 

Notice that at each time step, the numerical integration of the first integral over the fluid 
region in (8) is carried out by means of an isoparametric formulation which is used extensively 
in the finite element method, while the second integral can be evaluated analytically over 
the solid body B. The vorticity distribution over each distorted quadrilateral element in 
the physical space are actually performed over a square in the isoparametric space. Further 
details of the integration method can be found in Chen [5]. 

It is worthy of note that the evaluation of the integral (8) involves a problem associated 
with the non-uniqueness of the solution. The principle of vorticity conservation imposed 
by Wu [45] resolve such difficulty. For flows past single or multiple solid bodies, a more 
immediate improvement to the principle of vorticity conservation is provided by Chen [5]. 

Step 4 • The new internal velocity field can be established by solving the Poisson equations 
(11) with prescribed solid boundary conditions and outer boundary conditions of the velocity 
that have been determined from step 3. 

The final form of the discretized Poisson equations can be written 

{%:i ^ 

where u and v are two vectors of unknown interior nodal values. Also, A is a 11-banded 
matrix, while fi and {2 are vectors associated with the known forcing terms and boundary 
conditions. The resulting 11-banded matrix equations are then solved by a preconditioned 
biconjugate gradient routine [6]. This step completes the computational loop for each time 
level. 

One further important point to be noted in the integral approach is that the initial flow 
field can be determined by the same solution procedure described above (from step 2 to 
step 4). To be more precise, at time t = 0 + the unknown boundary vorticity values is 
determined by solving integral (9) with zero vorticity values everywhere away from the solid 
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boundary surface. This is based on the fact that the vorticity will concentrate on the body 
surface in the form of a vortex sheet immediately after the body impulsively started. This 
approach therefore reduces equation (9) to a boundary integral, indicate that the method 
can be viewed as the boundary element method utilized in the problem of the potential flow. 
Once the boundary vorticity values are obtained, the initial velocity field can be determined 
by solving equation (11) with the known velocity values at all points on the outer perimeter 
of the computational domain which have been calculated by the integral (8). In contrast to 
the special technique used by other methods, this integral approach enables the numerical 
code to generate the initial velocity field simply by the implementation of one cycle of a 
solution procedure rather than employing any additional treatments. 

3. RESULTS AND DISCUSSIONS 

In this section we apply the numerical scheme described above to simulate and control 
the unsteady flow around a rotating cylinder that undergoes a variety of steady and unsteady 
angular/rectilinear speed ratios at a Reynolds number of 200. The choice of this particular 
Reynolds number is not due to the limitation of the numerical algorithm, it is mainly for 
the purpose of comparing with the existing experiments of Coutanceau and Menard in the 
case of constant rotation for a rotating cylinder [8]. In this model, the rectilinear velocity is 
fixed as a constant value while the angular velocity is treated as a control variable. 

Although the choice of time-dependent rotation rates that may be used to control the 
rotating cylinder are unlimited, the computational results presented here are restricted to 
the following three types of rotation: 

1. Constant speed of rotation: a(t) = constant. 

2. Time-harmonic rotary oscillation: a(t) = Asin7rFt. 

3. Time-periodic rotation: a(t ) = A| sin ir(F/2)t\. 

Again, this choice was made because it allowed us to compare experimental and numerical 
results and it matched the control experiments of Modi [25, 24]. All variables are normalized 
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to the nondimensional forms in the formulation. In a non-rotating frame attached to the 
cylinder, the configurations for the different controls considered in the physical space are 
sketched in Figure 1, together with the corresponding time evolution of the angular velocity. 

Concerning the use of a time-harmonic rotary oscillation, it is common to define the 
motion of angular rotation 0(f) as 

0(t) = — 6 cos2n ft, (17) 

where 9 is the angular amplitude and / is the forcing frequency of the oscillation. Thus, the 
associated time-dependent speed ratio is given by 

. . afl(t) 

a(t) = — — — = AsmnFt, (18) 

where F = 2 af/U is the reduced forcing frequency and A = nF9 is the normalized maximum 
rotation rate of the forcing oscillation. In order to attain high lift and reduced drag, previous 
work for constant rotation rates [7] lead us to consider one special type of time-periodic 
rotation. That is, the cylinder under control is rotated in the counterclockwise direction 
about its axis with a time-periodic speed ratio given by 

a(t) = A \ sin n(F/2)t\. (19) 

Here, the reduced forcing frequency of this particular time-periodic rotation is F. This partic- 
ular type of rotation is expected to provide a substantial lift enhancement and drag reduction 
through a proper choice of both the angular amplitude (thus the normalized maximum ro- 
tation rate ^4) and forcing frequency (thus the reduced frequency F). This improvement can 
be demonstrated by comparing its respective force performances against the time-harmonic 
rotary oscillation. 

The major goal of this paper is to study the effect of rotation rate control upon the 
lift and drag on the cylinder surface. Hence, in the following discussion we concentrate on 
various issues concerning the development of temporal forces. In a viscous flow, it is well 
known that the total lift and drag forces are contributed by the pressure and skin friction due 
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to the viscous effects. An important consequence of using the velocity /vorticity formulation 
is that the forces can be directly evaluated from the known vorticity on the cylinder surface. 
Hence, for known vorticity values on the cylinder surface, the lift and drag coefficients can 
be calculated in the r-0 coordinates by 

C L {t) = C Lp {t) + C Lf (t) = J o cos0dd + Jfef o w (0r«* 6d0, (20) 

and 

Co(t) = Cr> p (t) + Cof{t) = — ^ sin 0 d9 — — jf U)(t) r sin9d9, (21) 

where the subscript T denotes quantities evaluated on the cylinder surface. The subscripts p 
and / represent the contribution from pressure and skin friction, respectively. In particular, 
we denote the positive values of Cl in the negative y-direction (as noted in Figure 1). 

To assess the accuracy of the numerical algorithm, computations were first performed 
over a wide range of constant speed ratios up to 3.25 at a Reynolds number of 200. Several 
particular speed ratio parameters were chosen to allow for the comparison against the ex- 
perimental work of Coutanceau and Menard [8]. Speed ratios greater than 2 are important 
in the study of the possibility of suppressing vortex shedding by an application of higher 
rotation rates. The details of the work using constant rotation rates were reported in [7], 
where several numerical solutions were compared and demonstrated to be in good agreement 
with experimental results. 

In the computations below, a fixed flow domain is used and its extent is essentially 
determined by the time-span under investigation. Namely, a larger computational domain 
will be needed for a longer time of observation. Here, a circle of radius r = 24 and r — 36 
are chosen for the time-span of 0 < t < 24 and 0 < t < 36, respectively. These long time 
histories were necessary to demonstrate the periodicity of the flow pattern and evolution of 
the forces. When a larger domain is considered, the mesh is increased in the radial direction 
in order to properly discretize the domain. A uniformly spaced M — 128 in circumferential 
direction is used for all computations, while N = 120 and N = 180 stretched grid lines 
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in the radial direction are used for the time-span under consideration. In the numerical 
calculations, small initial time steps axe taken in order to contain the inherent numerical 
error caused by resolving of the infinitesimal vorticity layer at t = 0 + . A detailed discussion 
of the algorithm and the accuracy of the initial flow field can be found in [7]. 

3.1 Constant Speed of Rotation 

Figures 2(a,b) show calculated instantaneous streamline plots for a constant value of 
speed ratio of a = 2.07 at Re = 200 and compare these plots with the experimental work 
of Coutanceau and Menard [8]. In the computation, the non-rotating reference frame is 
translating with the cylinder while the camera in the experiment is moving with the cylinder 
as well. Excellent agreement is obtained, despite the fact that a high velocity gradient is 
induced in the near wake due to the cylinder rotation. Although not shown here, one partic- 
ular interesting feature is the difference between the experimental work and our calculated 
observation regarding the conclusion of suppressing of vortex shedding at high speed ratios. 
In the computation using high speed ratios (at Re = 200, as reported in [7]), the calculated 
equi-vorticity contours seem to imply that vortex shedding continues to occur even at high 
rotation rates (a > 2.07). However, at these high a, the observed formation of the vortex 
street behind a rotating cylinder seems to contradict the experimental conclusion described 
in [8]. This difference is due to the fact that the experimental apparatus was such that only 
10 dimensionless time units of data could be collected and in part by the flow visualization 
techniques used in their experiments. On the other hand, it is important to note a recent 
investigation by Badr et al. [2] regarding the issue of suppressing of vortex shedding. Their 
tests were performed both experimentally and numerically at Reynolds numbers of Re = 10 3 
and Re = 10 4 . For a rotation rate of a = 3 at Re = 10 3 , they show that no other eddy is 
created after the shedding of two vortices. In addition, the temporal evolutions of the lift 
and drag coefficients imply that a steady state is indeed approached. 

Figure 3 shows plots of the time histories of lift, drag and lift/drag coefficients at various 
values of speed ratios (0 < a < 3.25) and for time in the interval 0 < t < 24. As seen in Fig- 
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ure 3(a), when the speed ratio is increased to 2.07, the lift increases timewise proportionally. 
However, as the speed ratio further increases, lift appears to initially decrease then increases 
gradually at later times. Not surprisingly, the maximum value of C L that can be achieved by 
rotation is also higher as the speed ratio grows. It is also observed that, at speed ratios lower 
than 2, the respective lift curves exhibit a well established periodic evolution. However, in 
the range of a > 2, it is not known whether the nature of this periodicity will continue 
if the time of investigation is expanded. Apparently, as can be seen from Figure 3(a), the 
cylinder rotation (as a boundary moving mechanism) does yield a substantial lift enhance- 
ment as indicated by the experimental work of Modi’s group for airfoil/rotating cylinders 
configurations. 

As illustrated by the drag curve in Figure 3(b), there is a substantial increase in drag 
when the speed ratio is increased. In all cases considered here, these drag curves seem to 
converge after a certain time and then oscillate under different amplitudes and frequencies 
thereafter. Detailed numerical results on the effect of the speed ratio to the resulting lift/drag 
curve are shown in Figure 3(c). In the range 0 < a < 2.07, the lift/drag performance appears 
to improve timewise (for 0 < t < 24) with an increase of a. If a comparison is made between 
a = 2.07 and a = 0.05, a noticeable improvement of the lift/drag performance is observed. 
Although a higher lift/ drag ratio is achieved by increasing the rotation rate in this range, the 
question arises whether any further increase of ot will result in a continued improvement of 
the lift/drag ratio. Intuitively, it is natural to expect a monotonical increase in the lift/drag 
ratio as o; increases to a = 3.25. However, this is not the case as a comparison is made 
between a = 3.25 and a = 2.07. In fact, the lift/drag curves illustrate a gradually decrease 
in performance over certain time interval when the speed ratio increases beyond 2. Moreover, 
this tendency toward lower lift/drag ratio becomes noticeable when a reaches the highest 
value (a = 3.25) considered here. Nevertheless, for all a considered here, a significant 
increase in the maximum value of Cl/Cd can be obtained by increasing a. However, it is 
found that it will reach its maximum value at a much later time for higher values of a. 
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3.2 Time-Periodic Rotation vs. Time-Harmonic Rotary Oscillation 

The previous results only applied to constant rotation rates. In this section we consider 
time- varying rotations. Although the maximum value of lift/drag ratio can be achieved at 
a certain optimal constant rotation rate described above, this does not immediately imply 
that the constant rotation will always yield a global maximum value of the lift/ drag ratio 
over a specific time interval. In fact, the general problem is to determine the optimal control 
input among all time dependent rotation rates. As a first step towards solving this problem, 
one needs to extend the computational calculations to the case of arbitrary time-dependent 
rotation rates. Since we shall not attempt to solve the necessary conditions corresponding to 
the optimal rotation rate problem and because the main goal of the paper is to gain insight 
into the possible form of an optimal controller, we shall concentrate on two periodic inputs 
and investigate the impact of each control on this problem. Hence, to keep the paper of 
reasonable length and the simulations simple we restrict our study to two periodic inputs. 

It is well known that when a cylinder oscillates in a uniform flow, the associated forcing 
oscillating frequency and amplitude can influence the vortex formulation and forces response 
substantially [42, 44]. It has been experimentally shown that at Re = 200, the natural 
Strouhal frequency of a non-rotating circular cylinder ( a = 0) is approximately F n = 0.185 
[43]. It is of important to study the behavior of fluctuating forces at imposed forcing fre- 
quencies which lie in the neighborhood of the natural frequency. The temporal evolutions of 
lift, drag and lift/drag are shown separately in Figures 4(a,b,c) for a time-periodic rotation 
a(t) = |sin0.25t| and a time-harmonic rotary oscillation a(t) = sin0.5f, respectively. In the 
case of time-periodic rotation, the cylinder under control is rotated in the counterclockwise 
direction about its axis with a time-periodic angular velocity. Notice that these two types 
of rotation are employed by the same forcing frequency (i.e. F = 0.16) which lies in the 
neighborhood of the natural frequency. The numerical results clearly confirm the expected 
benefit of this time-periodic rotation for both lift and drag forces, as shown in Figure 4. 

In comparing these two types of rotation, it should be noted that rotating in the same 
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direction causes the lift curve to be shifted upwards due to the nature of rotation, while the 
drag curve is shifted downwards. In terms of performance, this corresponds to an increase of 
the time-averaged lift force in the time-span of the investigation, while in the same time in- 
terval, a substantial reduction of the time-averaged drag as well. The resulting improvement 
of the lift/drag ratio is shown in Figure 4(c). 

To demonstrate the influence of this time-periodic rotation on the temporal development 
of these force coefficients, two additional values of the forcing frequency were tested. In 
Figure 5, all the time histories of lift, drag and lift/drag coefficients are shown for a forcing 
frequency (i.e. F = 0.08) which is 1/2 the frequency in Figure 4. As it can be seen from 
Figure 5(a), except in an initial stage, lift is increased for the time-periodic rotation («(*) = 

| sin 0.125<|) when compared to the time-harmonic rotary oscillation (<*(t) = sin0.25t). Also 
the drag is reduced as illustrated in Figure 5(b). However, the amount of drag reduction is 
not significant over the time-span of the investigation. As would be expected, the resulting 
lift/drag curves shown in Figure 5(c) exhibit almost the same behavior as those lift curves 
illustrated in Figure 5(a). 

Similar results are shown in Figure 6, where a higher rotation frequency (i.e. F = 0.32) is 
imparted to the cylinder. As before, the time-periodic rotation at this frequency enhances the 
lift performance when compared to the time-harmonic rotary oscillation. However, notice 
that the corresponding drag curves are oscillating about an “average value” in the time 
interval under consideration, for both types of rotation. The improvement in the lift/drag 
curve shown in Figure 6(c) is more noticeable when compared to the lift/drag curve in 
Figure 4(c). Notice that none of the frequencies considered in Figures 5 and 6 are in the 
neighborhood of the natural frequency. 

In a detailed examination of all comparisons between two types of rotations mentioned 
above, it is found that lift enhancement is not always timewise in the time interval under 
investigation. The exact time at which the time-periodic rotation outperform the time- 
harmonic rotary oscillation depends on the forcing frequency imparted to the cylinder. This 
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occurs earlier for a higher value of forcing frequency. However, the drag curves have a signif- 
icantly different character. The decrease or increase in drag depends on whether its forcing 
frequency lies in the neighborhood of natural frequency or not. If the forcing frequency is in 
the neighborhood of the natural frequency as the case in Figure 4(b), the drag curve exhibits 
a substantial reduction. On the other hand, for the forcing frequencies which are not in the 
neighborhood of the natural frequency as shown in Figures 5(b) and 6(b), the drag curves 
changes only slightly. In addition, for these forcing frequencies, the respective lift/drag 
curves are similar to the corresponding lift curves. This is due to the minor variations in the 
drag curves. 

3.3 Effects of Forcing Frequency and Angular Amplitude 

Because of the above observations, it is worthwhile to investigate the effects of vari- 
ous control parameters upon the force coefficients. A comparisons of force coefficients are 
shown in Figure 7, corresponding to all three forcing frequencies described above. It can be 
readily seen that these forcing frequencies have considerable influence on the amplitude and 
frequency of the oscillatory forces. 

In terms of performance, Figure 7(a) presents the evolution of lift coefficients and shows 
no clear advantage of changing the forcing frequency. However, as shown in Figures 7(a,b), 
the rotation ( a(t ) = | sin 0.25t|) which lies in the neighborhood of the nature frequency 
achieves a higher value of (C L ) max , and yields a slightly larger value of drag when compared 
to the other two cases (namely, a(t) = |sin0.125f| and a(t) = |sin0.5£|). The lift/drag 
curves shown in Figure 7(c) exhibit similar temporal evolutions. 

It is also interesting to study the effect of angular amplitude on the temporal evolution 
of forces while the forcing frequency is fixed as a constant. Figure 8 shows that resulting 
forces on the cylinder can differ significantly at different angular amplitudes for a(t) = 
A\ sin0.314f|. This rotation corresponds to a forcing Strouhal number of 0.2 which is in the 
neighborhood of the natural Strouhal number of 0.185. The angular amplitudes considered 
were A = 1.0,2.07 and 3.25. Apparently, as can be seen from these figures, a larger angular 
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amplitude definitely yields an incremental lift coefficient over the time-span of investigation 
(0 < t < 36). However, initially the drag increases with an increase of A, then after a certain 
time it oscillates with almost the same amplitude and frequency around an averaged value. 
Consequently, this leads to a substantial improvement in lift/drag with increasing A, as 
clearly shown in Figure 8(c). Nevertheless, the effect of angular amplitude is very noticeable 
when compared to the effect of the forcing frequency shown in Figure 7. 

3.4 Time- Averaged Value of Forces 

We also examined . the averaged values of the force coefficients over the time span of 
investigation as various control parameters are altered. The time-harmonic rotary oscillation 
(a(t) = sin nFt) described in §3.2 was considered with three different values of forcing 
frequency. The forces were averaged with respect to the time (0 < t < 24). Figure 9(a) 
shows the time-averaged values of lift, drag, lift/drag as the forcing frequency F is varied 
between 0.08 and 0.32. This figure shows that the forcing frequency has a considerable 
influence on the time-averaged values of the force coefficients. The local maximum values of 
time-averaged lift, drag and lift/drag ratios correspond to the forcing frequency which lies 
in the neighborhood of the natural frequency. This particular feature was also observed in 
the numerical results of Mo [23] where it was shown that the drag peak occurs at the forcing 
frequency equal to the natural frequency. 

As for the cases of time-periodic rotation, variations of time-averaged forces coefficients 
with respect to the forcing frequency are presented in Figure 9(b). As illustrated in the 
figure, a forcing frequency in this range (i.e. 0.08 < F < 0.32) has little effect on the 
time-averaged forces. Although the difference in time-averaged drag is minor, the forcing 
frequency which lies in the neighborhood of the natural frequency (F = 0.185) corresponds 
to a larger time-averaged drag and a smaller time-averaged lift. 

The effect of angular amplitude on the time-averaged values of lift, drag and lift/drag 
coefficients is shown in Figure 9(c), for a(t) = A|sin0.314f| averaged over 0 < t < 36. For 
A in the range, 1 < A < 3.25, all the time-averaged values are almost linearly proportional 
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to the angular amplitude. Significant increases in lift coefficients with increasing angular 
amplitude is particularly noticeable. This can be demonstrated by comparing the case A = 
3.25 with A = 1. It represents a 240% increment of lift performance. However, a slight 
increment in the drag coefficients with increasing angular amplitude is also observed. A 
moderate improvement of time-averaged lift/drag ratio is also seen. 

In the case of a constant speed of rotation, it is also worthwhile to study the effect of 
speed ratios on the time-averaged lift, drag and lift/drag coefficients. These results are 
shown in Figure 10(a) for ot in the range, 0 <[ ci ^ 3.25, and in the time interval 0 <! t ^ 24. 
It illustrates that the time-averaged lift is almost linearly proportional to the speed ratio, 
while the time-averaged drag remains as a constant value up to a = 2, then monotonically 
increases with speed ratio thereafter. Most importantly, the time-averaged lift/drag is not 
linearly proportional to the speed ratio. As shown in the figure, the highest value of the 
speed ratio a = 3.25 considered here is not the optimal constant rotation rate corresponding 
to the maximum value of time-averaged lift/drag. The maximum value occurs at a lower 
speed ratio, approximately a = 2.38, and it represents a substantial increase of 440% over 
the lower speed ratio a = 0.5. In Figure 10(b), the variation of the (total lift)/ (total drag) 
force ratio with respect to the speed ratio is shown for a in the range 0 < a < 3.25. Although 
a maximum is achieved at a value between c* = 2.0 and o; = 2.38, it should be noted that 
this optimal speed ratio is not necessarily the same optimal value as described in Figure 
10(a). 

The results presented in Figures 10(a,b) demonstrate an effective way of improving 
lift/drag performance by changing the rotation rate and illustrate the important of selecting 
a proper rotation rate in order to maximize performance. If one formulate this problem as 
an optimal control problem, then the time-averaged lift/drag ratio in Figure 10(a) represents 
the cost function and the goal is to find an optimal control aj , among a set of restricted con- 
trol parameters (constant values), that will provide the maximum value of the time-averaged 
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performance functional 
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where Tj is the final time. Similarly, the curve in Figure 10(b) represents the optimal control 
problem defined by maximizing the cost functional 
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In theory, these optimal control problems can be solved by applying necessary conditions 
for distributed parameter systems. However, computational methods for such necessary 
conditions are complex and not yet fully developed. The results presented here may be used 
to guide and test future computational schemes based on optimal control theory. 

3.5 Synchronizations of Cylinder and Wake 


The synchronization of cylinder and wake has long been known to be an important 
component of vortex-induced oscillations [33]. A detailed study of various types synchro- 
nization for a body oscillating transversely in a uniform stream can be found in Williamson 
and Roshko [44]. For the case of time-periodic rotation considered here, it is natural to ask 
whether such synchronization can occur and how well the numerical results can predict the 
occurrence of this important phenomenon. To the best of our knowledge, the current study 
is the first work to investigate synchronization under the particular form of time-periodic 
rotation described in previous sections. 

An examination of the responses in Figure 4, shows that the combined system of cylinder 
and wake will be “locked in” by an imposed forcing frequency. This synchronization of 
the cylinder and wake is due to the fact that the forcing frequency of rotation ( F = 0.16) 
lies in the neighborhood of the natural frequency ( F n = 0.185). Notice that in the case 
of time-periodic rotation shown in Figure 4, both lift and drag curves oscillate with the 
forcing frequency (corresponding to a time period of T = 12.5), clearly exhibiting a periodic 
response. However, in the case of time-harmonic rotary oscillation, the lift curve oscillates 
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with the same forcing frequency ( T = 12.5) while drag curve oscillates with the period of 
T/2. Consequently, the lift/drag ratios oscillate at the same frequency (T = 12.5) for both 
types of rotation. 

For the case of time-periodic rotation a(t) = A\ sin0.314<|, we extend our observation to 
a relatively longer time. For 0 < t < 36, an examination of these force curves for A = 1.0 
in Figure 8 exhibits a periodic response with a frequency ( F = 0.2) precisely equal to the 
input forcing frequency (i.e. T = 10). Although this periodic behavior is not established for 
A = 2.07 and 3.25, the corresponding curves are almost periodic in time. In order to confirm 
this periodicity, a sequence of instantaneous streamlines are shown in Figure 11. In Figure 
11, each plot is separated with an interval of one time period. These streamlines are plotted 
in a frame fixed with the undisturbed fluid. The periodicity of the flow is clearly noticeable. 
Two opposite-sign vortices are shed alternately on opposite sides of the cylinder at each cycle 
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of rotation. The vortex formation in the wake is similar to the case of a non-rotating cylinder 
(a = 0). However, the midline of the vortex street has been displaced slightly upwards due 
to the nature of rotation (in the counterclockwise direction). 

In order to identify the range of frequency for this fundamental synchronization, we im- 
pose a rotation rate a(£) = |sin0.283t| with a forcing frequency ( F = 0.18) which is the 
neighborhood of the natural frequency. The time histories of lift, drag and lift/drag coeffi- 
cients shown in Figures 12(a,b,c) clearly demonstrate the periodic behavior of the response. 
The corresponding streamline plots (at each instant) are presented in Figure 13. These results 
show that there exists a range of forcing frequencies in which fundamental synchronization 
will occur. However, the precise range of forcing frequencies leading to synchronization has 
not yet been determined (computational time is the limiting factor). 

In the case of time-harmonic rotary oscillation, the effects of the forcing frequency and 
amplitude on a cylinder wake have been investigated experimentally by Tokumaru and Di- 
motakis [40]. Several vortex formations were observed in the wake. Their experiments dealt 
with a range of amplitudes and frequencies at a Reynolds number of Re = 1.5 x 10 4 . By fixing 
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the reduced amplitude A in their experiments, four qualitatively different vortex shedding 
modes were identified when the forcing frequency was increased. Although their experiments 
were conducted at a Reynolds number higher than the current study (Re = 200), we con- 
ducted a similar investigation. Here, h particular value of the rotation rate a(f) = 2sin3.14f 
(corresponding to A = 2 and F = 1.0) was tested. Figure 14 shows the lift and drag histories 
up to t = 36. These curves are clearly periodic in nature. Under this forcing frequency, the 
lift curve oscillates with a period of T — 10, while the drag curve oscillates with a period 
of T = 5. The instantaneous streamlines plots are presented in Figure 15, and they show 
a time periodic flow pattern. Moreover, these results indicate that rotation may provide an 
effective control of the cylinder wake. 

4. CONCLUSIONS 

An algorithm for computing the viscous flow past a rotating cylinder is presented and 
applied to the problem of controlling cylinder forces by rotation. Several fundamental types of 
rotation were considered. Using time-periodic rotations leads to a considerable improvement 
in the force coefficients and was shown to be very effective, especially compared to time- 
harmonic rotary oscillations. These results are significant because they show a proper choice 
of the rotation rate can lead to improved flow fields. Very precise periodicity of the force 
for certain cases was established, and this periodic behavior has considerable impact on 
controlling the vortex formation in the cylinder wake. For the case of a constant speed of 
rotation, two optimal control problems were considered and solved computationally. 

These results demonstrate the feasibility of using boundary mechanisms for controlling 
unsteady flows, and consequently can be applied to enhance the performance. Using such 
mechanisms as a controller allows us to formulate a wide variety of optimal control problems 
for fluid flow systems. Modifications of existing numerical algorithms needed for such control 
problems depend on performance and design constraints. For example, one may need the 
maximum (or minimum) lift to drag ratio in order to sustain a particular maneuver of a 
supermaneuverable aircraft. Because of the complexity and importance of the relationship 
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between vortex motion and cylinder forces, the first step in control design may be to seek a 
specific type of rotation rate that will match all proposed goals. 

A precise understanding of time-dependent moving surfaces in boundary layer control 
may provide an effective way for lift enhancement and drag reduction. By treating the 
rotation rate as a control variable in this model, we will eventually be interested in finding 
the optimal control (i.e. a time history of the rotation rate) that maximizes (or minimum) 
the lift-to-drag ratio over a fixed time interval. Although here the optimal control problem 
associated with the constant rotation rate was solved by direct computations, it is still 
important to explore the possible implementation of a computational algorithm to calculate 
the optimal solution for the more general problems. The tools developed here can be used 
to investigate fundamental questions regarding control of separated flows by using various 
boundary control mechanisms. Future work needs to be done in the development of new 
computational algorithms for solving complex optimal flow control problems. 
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FIGURE 2: Instantaneous streamlines plots for Re = 200, a = 2.07 at (a) t = 5.0, (b) 
t = 9.0: computed (left); flow- visualization pictures (right). 
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FIGURE 3: Temporal evolution of the lift (a), drag (b) and lift/drag (c) coefficients at 
Re = 200 with various constant speed ratios (0.05 < a < 3.25). 
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FIGURE 4 : Comparison of temporal evolution of the lift (a), drag (b) and lift/drag (c) 
coefficients for a time-periodic rotation a(t) = | sin 0.25i| (T = 12.5) with a time-harmonic 
rotary oscillation a(t) = sin0.5< (T = 12.5) at Re — 200 for 0 < t < 24. 
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FIGURE 6: Comparison of temporal evolution of the lift (a), drag (b) and lift/drag (c) co- 
efficients for a time-periodic rotation a(t) = | sin 0.5t| with a time-harmonic rotary oscillation 
a(t) = sinf at Re = 200 for 0 < t < 24. 












FIGURE 8: Temporal evolution of the lift (a), drag (b) and lift/drag (c) coefficients 
for a time-periodic rotation a(t) = A|sin0.314t| at Re = 200 with various amplitudes of 
A = 1.0,2.07 and 3.25 for 0 <t< 36. 
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FIGURE 10: Effect of the speed ratio on time-averaged lift, drag and lift/drag coefficients 
(a) and on total lift/total drag force ratio (b) for 0 < a < 3.25. 
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FIGURE 11: Instantaneous streamlines plots for Re = 200, a(t) = |sin0.314f| (F = 0.2), 
viewed from a frame fixed with the undisturbed fluid, (a) t = 16, (b) t = 26, (c) t = 36. 
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FIGURE 13: Instantaneous streamlines plots for Re = 200, a(t) = | sin0.283t| ( F = 0.18), 
viewed from a frame fixed with the undisturbed fluid, (a) t = 14, (b) t = 25, (c) t = 36. 
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(a) 



(b) 



FIGURE 14: Temporal evolution of the lift (a) and drag (b) coefficients for a time-harmonic 
rotary oscillation a(t ) = 2sin3.14f ( F = 1.0 and A = 2) at Re = 200 for 0 < t < 36. 
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FIGURE 15: Instantaneous streamlines plots for Re = 200, ce(t) = 2sin3.14f (F = 1.0 
and A = 2), viewed from a frame fixed with the undisturbed fluid, (a) t = 16, (b) t = 26, 
(c) t = 36. 
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